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Abstract.  In  this  paper  we  report  our  recent  results  on  the  design  of  materials  with 
controlled  properties  by  the  application  of  computational  chemistry  methods,  new 
algorithms  and  scalable  software. 

1  Introduction  ;  62_"!17 

Developments  of  new'  materials  (such  as  new  high  energy  density  mate 

oligomeric  silsesquioxanes)  that  are  highly  resistant  to  extreme  environ 

desirable  coatings  for  rocket  engines;  nonlinear  optical  materials;  a 

important  in  many  applications  of  interest  to  the  Air  Force.  However, 

complex  material  systems  with  controlled  properties  is  challenging.  Th<  ,  —  wwv'uiauUUS 

we  describe  in  this  paper  would  not  have  been  possible  without  scalable  quantum  chemistry 

codes  such  as  GAMES  S  (General  Atomic  and  Molecular  Electronic  Structure  System), 

developed  under  the  auspices  of  a  DoD  CHSSI  (Computational  High  Performance  Scalable 

Software  Initiative)  grant,  Gaussian,  and  molecular  dynamics  codes,  as  well  as  the  availability  of 

the  DoD  Computation  Centers  through  our  computational  challenge  project 

2  Results  and  Discussion 

2.1  High-Energy  Density  Materials  (HEDM) 


2.1.1  Synthesis  of  new  polynitrogen  compounds. 

The  identification,  development,  and  formulation  of  new  energetic  materials  for  advanced 
rocket  propulsion  applications  is  an  area  of  long  standing  interest  to  the  Air  Force.  The 
performance  limits  of  currently  used  propellants  have  been  reached,  so  new  energetic 
compounds  are  required  to  significantly  improve  the  ability  of  the  warfighter  to  access  and 
control  space. 


Polynitrogen  species  such  as  the  recently  discovered  Ns"  cation  are  of  interest  as  potential 
energetic  ingredients  in  new  propellant  formulations.  The  recent  successful  synthesis  of  Ns+  in 
macroscopic  quantities  has  prompted  the  search  for  additional  polynitrogen  compounds. 
Computational  chemistry  plays  a  central  role  in  determining  the  stabilities,  potential  synthetic 
pathways,  and  key  spectroscopic  "fingerprints"  of  new  polynitrogen  species. 

The  structures,  stabilities,  vibrational  frequencies,  and  infrared  intensities  of  several  potential 
synthetic  precursors  to  new  polynitrogen  species  have  been  computed  using  ab  initio  electronic 
structure  theory,  at  the  second  order  pertuibation  theory  level  (MP2,  also  known  as  MBPT(2)), 
using  the  6-31G(d)  valence  double-zeta  polarized  basis  set.  Shown  in  the  figure  below  is  the 
predicted  structure  of  triphenylmethyldiazonium  cation,  also  known  as  trityldiazonium,  which  is 
a  possible  precursor  to  new'  polynitrogen  compounds  such  as  pentazole,  a  5-membered  ring 
system.  The  calculated  structure  shows  that  this  cation  is  unstable  with  respect  to  dissociation  of 
N2.  Therefore,  these  calculations  predict  that  this  cation  is  not  a  viable  polynitrogen  precursor. 


Figure  1.  MP2/6-31G(d)  optimized  structure  of  the  trityldiazonium  cation,  [CmHisNz]4-. 

Although  this  is  a  negative  result  in  the  sense  that  it  indicates  that  trityldiazonium  is  not  a 
stable  precursor,  it  is  nonetheless  a  highly  useful  result  in  multiple  ways.  First,  it  saves 
significant  time  and  effort  by  eliminating  from  consideration  for  subsequent  attempts  at  synthesis 
a  compound  which  is  not  likely  to  be  stable.  Furthermore,  it  suggests  ways  in  which  the 
trityldiazonium  cation  may  be  chemically  modified  in  order  to  overcome  its  instability  (e.g.,  by 
judicious  placement  of  electron-withdrawing  groups  on  the  trityl  moeity.)  Future  work  in  this 
area  will  examine  such  derivatives  of  the  trityldiazonium  cation,  as  well  as  other  classes  of 
promising  polynitrogen  precursors. 


2.1.2  High-energy  density  monopropellants 

A  recent  major  focus  centers  around  a  series  of  proposed  monopropellants[l],  such  as: 


Isp,  which  is  a  measure  of  the  energy  released  by  the  reaction  of  the  fuel  divided  by  the  mass  of 
the  products,  is  of  primary  interest: 

Isp  oc  [AH/mf2  (1) 

Using  a  combination  of  isodesmic  reactions  and  the  G2  model,  the  heat  of  formation  for  I  is 
predicted  to  be  456.8  kcal/mol.  This  translates  to  an  1^  of  329  sec,  compared  with  230  sec  for 
hydrazine,  thus  being  a  very  promising  fuel.  However,  it  is  important  to  consider  the  stability  of 
such  high-energy  species  to  various  reactions,  before  asserting  their  viability  as  fuels.  One 
possible  reaction  of  I  is  the  loss  of  molecular  nitrogen  to  form  the  smaller  species  shown  below: 


Using-  our  highly  scalable  second  order  perturbation  theory  code  with  a  large  basis  set,. we  have 
found  that  this  decomposition  process  is  exothermic  by  60  kcal/mol.  It  may  be  that  the  barrier  for 
decomposition  is  large,  in  which  case  the  species  may  still  be  viable.  Dr.  Jeff  Bottaro  (SRI) 
suggested  that  I  may  be  unstable,  because  it  could  open  at  either  or  both  ends  to  an  azide.  He 
suggested  instead  that  we  consider  an  isomer  of  I,  in  which  the  CNO2  group  on  the  left  in  the 


figure  is  shifted  one  position  to  the  left,  and  the  CNO2  group  on  the  right  is  shifted  one  position 
to  the  right,  giving  H  (shown  above),  since  this  structure  is  less  likely  to  open  to  an  azide  or  a 
diazide.  We  therefore  initiated  MP2  geometiy  optimizations  for  both  II  and  the  diazide  III.  We 
find  that  II  is  indeed  lower  in  energy  than  I,  by  about  15  kcal/mol,  but  that  the  diazide  HI,  in 
which  both  ends  of  I  are  opened,  is  much  higher  in  energy  than  either  I  or  EL  So,  it  appears  that 
II  should  be  studied  further  as  a  possible  HEDM,  but  that  I  may  still  be  viable  as  well.  We  are 
now  investigating  feasible  decomposition  products  of  both  I  and  H,  as  well  as  the  heat  of 
formation  of  II. 

Recently  the  first  new  all-nitrogen  compound  in  nearly  a  century,  N5+[2],  was  synthesized.  The 
counterion  for  this  species  was  AsFg',  the  result  being  a  rather  unstable  solid.  In  order  to 
understand  the  potential  energy  surface  for  this  species,  we  performed  MP2  calculations  on  the 
various  stationary  points.  The  first  step  in  this  process  was  to  separate  the  two  ions  by  a  large 
distance  and  optimize  the  geometry.  The  result  was  an  ion  pair  fNs+][AsF<f][  complex  that  is,  not 
surprisingly,  98  kcal/mol  below,  the  separated  ions.  What  was  somewhat  surprising  is  that  this 
complex  is  47  kcal/mol  below  the  neutral  species  that  is  formed  when  one  F'  ion  is  transferred  to 
Nj+  to  form  neutral  AsF5  +  FN5.  The  former  is  a  well-known  species,  but  the  latter  has  not 
previously  been  reported.  FN5  has  several  isomers,  whose  relative  energies  and  the  barriers 
separating  them  were  determined-  using  large  basis  set  CCSD(T)  calculations.  The  lowest  lying 
isomers  were  found  to  have  small  energy  barriers  (<  10  kcal/mol)  separating  them  from  the  much 
more  stable  decomposition  products  FN3  +  N2.  The  exothermicity  of  this  decomposition  is  more 
than  50  kcal/mol,  so  the  products  AsFs  +  FN3  +  N2  are  slightly  lower  in  energy  than  the  original 
ion  pair  [Nj'jfAsFg'].  Since  FN3  is  itself  an  energetic  species  with  a  relatively  small  barrier  to 
further  decomposition,  it  is  therefore  not  surprising  that  the  original  species  made  by  Christe  et 
al.  is  unstable.  In  order  to  study  this  system  further,  we  have  employed  the  dynamic  reaction 
path  (DRP)  method  in  GAMESS  to  analyze  the  effect  of  putting  photons  into  vibrational  modes 
that  might  lead  to  decomposition.  Interestingly,  this  molecule  appears  to  be  very  RRKM-like, 
since  the  energy  dumped  into  specific  modes  rapidly  gets  distributed  to  many  other,  apparently 
strongly  coupled,  modes.  So,  one  must  put  much  more  energy  into  these  modes  than  the  barrier 
height  in  order  to  induce  dissociation. 


2.1.3  Quantum  molecular  dynamics  (MD)  simulations  of  low  temperature  HEDM 

Traditional  classical  MD  algorithms  can  be  extended  to  incorporate  equilibrium  quantum 
mechanical  effects  through  the  use  of  discretized  Feynman  path  integrals[3].  The  resulting 
algorithm,  path  integral  MD  (P1MD),  is  highly  amenable  to  parallelization  because  each  physical 
particle  becomes  rigorously  mapped  onto  a  collection  of  classical-like  quasiparticles  at  different 
values  of  imaginary  time,  each  connected  to  its  neighbor  by  a  harmonic-like  force.  The 
discretized  PIMD  algorithm  can  be  load-balanced  across  computational  nodes  because  the 
quasiparticles  with  the  same  imaginary  time  index  experience  identical  forces.  The  inter-node 
communications  are  minimal  since  there  are  only  nearest-neighbor  harmonic  forces  to  compute. 
The  PIMD  algorithm,  no  matter  how  complex  the  physical  system,  can  always  be  made  to 
achieve  excellent  Scalability  and  thus  high  performance.  Moreover;  classical  MD  techniques 
such  as  constant  temperature  and  pressure  algorithms  can  be  easily  incorporated  into  PIMD.  The 
algorithm  is  even  amenable  to  a  two-tier  level  of  parallelism,  the  first  being  over  the  imaginary 


time  slices  (quasiparticles),  the  second  over  the  force  loop  for  the  interparticle  interactions  in 
large  (thousand  or  even  million)  particle  simulations. 

In  an  important  theoretical  advance[4-7],  PIMD  methods  have  been  extended  to  include  quantum 
dynamical  effects  (i.e.,  not  just  equilibrium  properties).  This  computational  approach  is  called 
"Centroid  Molecular  Dynamics"  (CMD).  It  incorporates  the  dominant  quantum  dynamical 
effects  of  a  many-body  system  into,  a  classical-like  MD  framework,  thus  significantly  extending 
the  range  and  applicability  of  MD  methods.  The  basic  result  of  CMD[5-8]  is  that  the.  dynamical 
correlations  of  quantum  particles  in  a  general  many-body  system  can  be  accurately  computed  by 
running  classical-like  trajectories  on  an  effective  potential  which  includes  the  effects  of  quantum 
zero  point  energy  and  tunneling.  The  task  of  integrating  the  CMD  equations  is  not  trivial  since 
the  effective  potential  is  a  quantum  potential  of  mean  force,  requiring  an  "on  the  fly"  dynamical 
averaging  at  each  timestep.  A  number  of  powerful  algorithms  for  solving  the  CMD  equation 
have  been  developed[8].  The  best  CMD  algorithm  so  far  is  the  “hyperparallel”  CMD 
approach[8]  using  an  adiabatic  algorithm  in  the  first  tier  of  parallelism  to  determine  the  effective 
potential;  a  second  tier  of  parallelism  can  be  introduced  in  which,  for  each  of  the  averaging 
slaves,  the  MD  steps  are  sent  to  additional  nodes  for  each  discretized  path  integral  imaginary 
time  slice.  A  third  tier  of  parallelism  for  very  large  systems  can  be  utilized  for  each  imaginary 
time  slice,  in  which  the  classical  force  loop  is  split  over  several  nodes.  The  scaling  of  the 
hyperparallel  CMD  code  on  the  IBM  SP  and  SGI  02000,  as  well  as  a  Pentium  Linux  cluster,  is 
excellent  (Fig.  2).  The  parallelization  of  the  long  rage  forces  (Ewald  summation)  was  also 
achieved. 


Cade  Safing 


Figure  2.  Scaling  of  CMD 

The  PIMD  and  CMD  algorithms  are  of  general  interest  in  all  areas  of  condensed  matter  computer 
simulation.  However,  they  are  particularly  useful  for  the  HEDM  program  in  the  computational 
study  of,  e.g.,  solid  hydrogen,  helium,  or  nitrogen  matrices  which  contain  isolated  reactive 
species  such  as  lithium,  boron,  aluminum,  hydrogen  atoms,  organic  radicals,  etc;  Such  matrices, 
are  high  priority  targets  as  possible  HEDMs  for  Air  Force  space  propulsion  purposes.  At  the 
temperatures  appropriate  to. the  condensed  phases  of  hydrogen  or  helium,  the  host  matrix 
molecules  will  exhibit  large  quantum- mechanical  effects,  but  the  treatment  of  these  effects  by  a 
direct  attack  on  the  time-dependent  Schrodinger  equation  is  impossible.  The  important  issue  of 
the  HEDM  is-its  stability.  It  must  be  able  to  trap  the  energetic  impurities  for  some  period  of  time 
and  thus  impede  their  diffusion  and  eventual  recombination.  Both  the  structural  and  dynamical 
aspects  of  the  composite  condensed  phase  HEDM  systems  are  important  to  understand.  For  the 
equilibrium  structural  studies  and  thermodynamic  studies,  PIMD  methods  are  implemented.  For 
the  dynamical  studies,  the  CMD  approach  must  be  employed  to  directly  simulate  the  diffusive 
and  recombination  steps  of  the  guest  atoms  in  the  solid  host.  The  long-term  goal  of  this  research 


is  to  develop  general  computational  methods  to  rapidly  and  efficiently  characterize  proposed 
high  energy  density  materials  and  to  draw  some  conclusions  about  the  utilizability  of  specific 
HEDM  systems. 

The  structure  and  stability  of  atomic  impurities  trapped  in  solid  para-hydrogen  have  been 
studied  by  employing  large  scale  PIMD  and  CMD  simulations  at  4  K  and  zero  external  pressure. 
Doped  systems  were  prepared  by  substituting  impurity  atoms  for  hydrogen  molecules.  The 
structural  and  thermodynamic  properties  of  the  resulting  HEDM  were  then  calculated.  In  the  case 
of  atomic  aluminum  impurities,  the  impurity  atom  interacts  anisotropically  with  the  matrix  para- 
hydrogen  molecules  because  of  its  singly  filled  3p  orbital.  To  assess  the  effect  of  the  electronic 
anisotropy  in  the  potential,  a  comparison  has  been  made  between  simulations  in  which  the 
orientation-dependent  (anisotropic)  and  orientation-averaged  (spherical)  AI-H2  potentials  were 
used.  In  collaboration  with  Millard  Alexander[9],  three  matrices  were  investigated:  (a)  one  with 
a  single  A1  atom  site-substituted  for  a  p-Hj  molecule,  (b)  one  with  a  similar  site-substituted 
matrix  having  a  nearest-neighbor  vacancy  (defect),  (c)  and  one  with  a  periodic  "super  cell" 
simulation  box  having  two  hydrogen  bulk  slab  regions  in  order  to  study  a  A1  impurity  near  the p- 
H2  surface, 

It  was  found  that  in  the  presence  of  an  adjacent  defect,  the  aluminum  moves  toward  the  defect, 
and  in  the  case  where  the  orientation  dependent  Al-p-Hfe  potential  is  used,  the  A1  occupies  the 
position  halfway  between  its  original  substitution  site  and  the  vacancy  site.  The  presence  of  the 
vacancy  helps  to  stabilize  the  A1  system  when  the  orientation  dependent  potential  is  used, 
whereas  it  was  found  to  have  little  effect  in  the  case  of  B  impurity[10].  Hence  the  aluminum 
atom  seems  to  prefer  occupying  a  double  substitutional  site.  The  effect  that  this  has  on  the 
dynamical  stability  of  the  A1  doped  system  will  be  studied  in  future  research. 

2.2  POSS  Compounds 

There  is  great  interest  in  POSS  compounds  because  of  their  resistance  to  extreme  environments 
and  the  ease  with  which  they  are  synthesized.  Since  almost  nothing  is  known  experimentally 
about  the  mechanism  by  which  they  are  formed,  we  have  embarked  on  a  long  term  project  to 
determine  the  possible  competing  mechanisms  as  a  function  of  catalyst,  solvent,  and  substituents 
(e.g.,  in  {11]),  The  synthesis  begins  with  hydrolysis  of  RSiX3  to  RSi(OH)3,  followed  by 
condensation  to  the  siloxane  RSi(OH)rO-Si(OH)2R  Subsequent  condensations  lead  to  the  3D 
cage  compounds.  We  have  shown  that  the  initial  hydrolysis  and  condensation  steps  all  have  very 
high-energy  barriers  in  the  gas  phase,  but  are  being  reduced  to  nearly  zero  by  the  presence  of  one 
water  molecule  added  to  represent  the  solvent.  Thus,  all  steps  leading -to  the  initial  disiloxanes 
and  to  the  ring  compounds  D3  and  D4  occur  with  net  energy  requirements  of  less  than  10 
kcal/mol.  In  this  paper  we  describe  current  efforts  to  (a)  determine  substituent  effects  on  these 
barriers,  (b)  determine  the  effects  of  adding  additional  water  molecules,  and  (c)  to  compare  the 
properties  of  the  3D  POSS  compounds  with  their  Ti  analogs.  The  substituents  studied  include 
R=H,  methyl,  t-butyl  and  phenyl,  and  X=CI,  OCH3.  There  is  experimental  interest  in 
synthesizing  incompletely  condensed  POSS.  Because  the  reaction  is  very  fast,  it  always  goes  to 
completion.  The  interest  in  substituent  effects  lies  in  attempting  to  slow  down  the  reaction, 
possibly  using  bulky  R  groups.  There  is  industrial  interest  in  the  possibility  that  Ti  POSS  may  be 
another  class  of  materials  with  desirable  properties.  Hammes-Schiffer  et  al.[12]  have  been 


developing  methodologies  for  the  simulation  of  hydrogen  transfer  reactions,  and  most  recently  a 
new  molecular  orbital  method  for  the  simultaneous  calculation  of  electronic  and  hydrogen 
vibrational  wavefunctions,  now  incorporated  into  GAMESS.  This  method  may  be  used  to  obtain 
minimum  energy  reaction  paths  and  direct  dynamics  trajectories  of  hydrogen  transfer  reactions, 
with  the  advantage  that  nuclear  quantum  effects  such  as  zero  point  energy  and  hydrogen 
tunneling  are  incorporated  during  the  generation  of  the  reaction  paths  and  trajectories,  rather  than 
subsequently  included  as  corrections.  This  approach  will  be  combined  with  a  mixed 
quantum/classical  surface  hopping  method  to  study  the  quantum  dynamics  of  hydrogen  transfer 
reactions  with  an  ab  initio  potential  energy  surface  obtained  “on  the  fly,”  and  used  to  investigate 
hydrogen  transfer  reactions  in  the  hydrolysis  and  condensation  steps  required  for  the  synthesis  of 
POSS,  in  order  to  determine  the  influence  of  different  trihalosilane  reactants  on  the  rates  and 
yields  of  these  steps  and  aid  in  the  efficient  synthesis  of  POSS. 

23  Nonlinear  Absorbing  Materials 

We  have  previously  shown  that  density,  functional  theory  and  the  time-dependent  density 
functional  response  theory  (TDDFT)  (B3LYP)  accurately  predict  the  structures  and  spectra  of 
porphin  and  the  larger  tetraphenylporphyrin  and  its  (3-octahalogenated  derivatives,[13,14]  thus 
enabling  their  assessment  as  reverse  saturable  absorbing  materials.  For  porphin  (PH,),  the 
photoinduced  tautomerization  process  has  also  been  investigated[15]  in  order  to  discern  between 
previously  proposed  discordant  mechanisms.  Recently  we  carried  out  calculations!  16]  for 
ietrabenzoporphyrin  (ZnTBP)  and  phthalocyanine  (ZnPc)  (Figure  3). 


Figure  3.  Zinc  phthalocyanine 

Our  calculations  support  the  multiple  bands  in  the  B  and  higher  energy  regions  of  the  ZnPc 
and  ZnTBP  spectra  that  have  been  determined  from  experiment.  Despite  considerable 
experimental  and  theoretical  /work  on  ZnPc  and  ZnTBP,  their  ground  state  electronic  structures 
and  UV  spectra  are  not  well  understood.  Previous  theoretical  studies  provided  only  a  qualitative 
understanding  in  the  low-energy  region,  and  none  for  the  higher  energy  bands.  We  have 
elucidated  significant  shifts  in  Gouterman’s  four  orbitals  of  ZnP  upon  tetraazasubstitution  and 
tetrabenzoannulation.[16]  The  benzo  groups  destabilize  the  laIn  orbital  in  ZnP,  but  they  do  not 
show  a  significant  deviation  from  the  four-orbital  model  for  the  interpretation  of  the  Q  and  B 
bands  of  ZnTBP.  However,  tetraazasubstitutions  considerably  lower  the  48^  orbital  (HOMO-1) 
in  ZnTBP.  Also,  it  is  shown  that  the  near  degeneracy  of  the  HOMO  and  HOMO-1  in  ZnP  that 
provides  the  basis  for  Gouterman’s  four  orbital  model  does  not  hold  for  ZnPc.  In  ZnPc,  the  near 
degeneracy  of  the  HOMO-1  (4aju)  with  other  occupied  orbitals  gives  rise  to  a  complex  structure 


of  the  higher  energy  regions  of  the  spectrum  that  were  not  reproduced  by  simple  MO  methods. 
The  4aju  MO  has  significant  %  contributions  from  the  benzo  rings  in  ZnPc  but  not  in  ZnTBP.  The 
Q  and  B  bands  are  predicted  to  appear  at  2.18  and  3.28  eV  for  ZnTBP,  respectively,  in  excellent 
agreement  with  the  maxima  located  at  2.02  and  3.18  eV  as  established  from  the  absorption 
spectrum  in  a  supersonic  jet  expansion,  with  similar  values  (2.02  &  3.10  eV)  obtained  from  the 
absorption  spectrum  in  an  Ar  matrix.  Tetraazasubstitutions  are  predicted  to  redshift  the  Q  band 
of  ZnTBP  by  about  0.1  eV.  The  computed  excitation  energies  in  near  quantitative  agreement 
with  experiment  clearly  provide  a  strong  basis  for  the  interpretation  of  the  electronic  spectra  of 
such  large  molecular  systems.  TDDFT  calculations  are  also  being  carried  out  for  candidate  two- 
photon  absorption  (TP A)  materials, [17]  to  be  compared  to  those  of  Albota  et  al.[18] 

2.4  Liquid  Crystals 

Recent  work  is  in  progress  for  studying  a  liquid  crystalline  droplet  to  model  the  behavior  of  a 
liquid  crystal  in  the  bulk,  with  the  details  reported  elsewhere. [19]  Indeed,  atomic  level 
simulations  are  proven  important  in  understanding  the  structure-property  relations  of  materials. 
Due  to  the  size  of  the  LC  droplet,  very  large-scale  classical  MD  simulations  of  such  systems  are 
required,  that  may,  become  computationally  prohibitive.  The  fast  multipole  method  (FMM)[20], 
which  uses  a  multiscale  hierarchy  of  partitions  of  the  volume  and  a  divide-and-conquer  strategy 
to  compute  the  power  series,  allows  all  the  forces  to  be  computed  to  any  specified  accuracy  in 
0(N)  operations.  FMM3D  (ah  implementation  and  improvement  to  FMM  in  3-D),  which 
contains  a  variety  of  schemes  for  computing  multipole  translations,  has  been  implemented  into 
the  MD  program  that  is  being  developed  in  our  group.[21] 
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